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COMPARISON OF POST-STORM NON- ADIABATIC RECOVERY 
OF TRAPPED PROTONS WITH RADIAL DIFFUSION 

Finn So'raas* 

ABSTRACT 

During magnetic storms it has been observed that the geomagnetically 
trapped protons undergo rapid non-adiabatic changes followed by slow non- 
adiabatic recovery approaching the pre-storm values. The slow non-adiabatic 
recovery can be accounted for in a semiquantitative way by solving a time 
dependent Fokker-PIanck equation with radial transport and loss terms describ- 
ing coulomb energy degradation and charge exchange. The equation is solved 
numerically in a region of space where we have measurements of the 100 keV to 
1700 keV protons mirroring at the equator. The transport term is assumed to 
have the form D = k/J" L p where is the magnetic moment of the proton and L 
the Mcllwain shell parameter. The value of D which gives the best fit to the 
data is found. Due to the limited amount of data used in this study the m and L 
dependence of the diffusion coefficient are not determined very accurately. The 
resulting values of D was found to be larger than the value evaluated by Nakada 
and Mead. The e -folding time for the intensities of the higher energies to 
recover at L = 3.5 is of the order of a year. 


*NAS-NASA Postdoctoral Resident Research Associate 
Present address: Department of Physics, University of Bergen, Bergen, Norway. 
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COMPARISON OF POST-STORM NON -ADIABATIC RECOVERY 
OF TRAPPED PROTONS WITH RADIAL DIFFUSION 

INTRODUCTION 

A short time alter the discovery of the geomagnetically trapped particles it 
was suggested by Kellogg (1959) that the spatial and temporal behavior of these 
particles might be influenced by processes that produced departure from first 
order adiabatic motion, resulting in diffusion across magnetic shells, conserving 
the first two adiabatic invariants and violating the flux invariant. Such a process 
has been effectively employed by Nakada and Mead (1965) in an attempt to arrive 
at an equilibrium outer zone proton distribution. 

An apparent radially inward motion of MeV electrons has been observed by 
Frank (1965) near the equatorial plane and at high latitudes by McDiarmid and 
Burrows (1966) and by Craven (1966). One possible interpretation of these data 
invokes trans-L diffusion. Newkirk and Walt (1968) have evaluated the radial 
diffusion coefficient for electrons at low L-values and Kavanagh (1968) for elec- 
trons in the region four to seven earth radii. 

Sbraas and Davis (1968) have made a study of the temporal behavior of the 
100 keV to 1700 keV trapped protons in the range 2 to 5.5 earth- radii for the 
period 29 of January through 29 of June 1965. They found that the particle in- 
tensities exhibit both adiabatic changes which varied directly with the magnetic 
field (Dst) and non-adiabatic changes which do not track the field. They were 
able to transform the proton intensities measured in the time variable field to a 
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reference dipole field using a transformation which conserved the three adiabatic 
invariants of motion. The transformed intensities then supposedly show time 
variations due only to non-adiabatic processes in which at least one invariant 
was violated. 

Their results show two types of non-adiabatic variations. First during 
magnetic storms the protons underwent rapid changes. Protons having energies 
less than about 200 keV were enhanced while higher energy protons were de- 
pleted. During the 18 April 1965 storm, which was the largest for this period, 
the enhancement was as much as a factor of four and the depletion as much as a 
factor of 10. The second type was a /low post-storm recovery of both high and 
low energy protons toward their prestorm values. The recovery times were 
observed to decrease with increasing radial distance. They also observed that 
the proton energies before the large april 18 storm exhibited the L~ 3 dependence 
characteristic of cross- L diffusive equilibrium and did not show this after the 
storm, though they slowly returned toward this dependence. These last two re- 
sults coupled with the results of Nakada and Mead (1965) suggest that the recover 
phase might be governed by radial diffusion of the protons. 

In the present analysis the temporal behavior of the 100 keV to 1700 keV 
protons between the April 18 and June 15, 1968 geomagnetic storms will be 
compared with a calculated time-dependent radial diffusion model and used to 
evaluate the radial diffusion coefficient. 
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DIFFUSION CALCULATION 

The equation used to describe the motion across L- shells of protons mir- 
roring at the equator is a one-dimensional Fokker- Planck equation with radial 
transport terms and loss terms describing coulomb energy degradation and 
charge-exchange. 


8n S , . 1 § 2 , . 8 //Afx\ \ n 

It = "57 ( D , n ) + ( D 2 n ) " V \\Ki/ n ) ' r 


( 1 ) 


where n (r, ,u, J, t) dr du dj is the number of particles between equatorial radius 
r and ( r + dr), having magnetic movement between p- and (p + 4 a ) ) and mtegral 
invariant J to (J +dj) at the time t . The coefficient is the mean radial 
displacement per unit time ^Ar/At^ and D 2 is the mean square radial displace- 
ment per unit time ^(Ar) 2 /At^. The angular brackets denote time averages. 

The third term on the right hand side of the equation describs coulomb energy 
loss and the fourth term charge exchange loss, r is the e -folding lifetime for 
charge exchange. In the case of coulomb energy loss the higher energy particles 
are a source for low energy particles whereas charge exchange removes parti- 
cles catastrophically. Any pitch angle scattering of the protons are neglected. 
The equation is the same as used by Nakada and Mead (1965) and similar to those 
used by Davis and Chang (1962), Tverskoy (3 564), and Newkirk and Walt (1968) 
in studies of trapped particles. 
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Falthamm ar (1966) has shown that D x and D 2 are related by 


D 


1 




( 2 ) 


where 

D = D 2 /2 . 

As a consequence the Fokker- Planck equation (1) can be rewritten as a diffusion 
equation of the form 



Falthammar (1967) has found a general expression for the diffusion coefficient D 
which depends on the power spectrum of the electromagnetic disturbances. He 
points out that if the magnetic disturbances have a power spectrum falling off 
with the frequency raised to the negative exponent of s then the diffusion coeffi- 
cient is of the form 


b u r 6+2s • (J . 2-9 (4) 

where the subscript M on D is used to indicate that this form of the diffusion 
coefficient is due to magnetic disturbances and not to time dependent electric 
potential fields. 

The diffusion coefficient is here assumed to have the form 

D = k • • r p (5) 
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By choosing different values of k, m and p Equation (3) can be solved for differ- 
ent values of the diffusion coefficient. 

For the coulomb energy loss and charge- exchange lifetimes entering into 
Equation (3) the same expressions as used by Nakada and Mead (1965) are used. 
They are 


aA , ri) 

fat ) - - 3. 55 * 10 6 P y ^ J MeV/ Gauss • day 

with r in earth radii, p in electrons/cm 3 and p the magnetic moment of the 
protons in MeV/Gauss. The electron density is given by 



p ~ (8000/r 4 +50) electrons/' cm 3 

the charge exchange term is 1/r = o-p Q v , where p 0 is the neutral hydrogen 
density, given by 


P Q = 7.35 * 10 3 /r 5 


atoms/cm 3 


For the charge excliange cross section, cr, values *?*ven by Allison (1958) were 
used v is the proton velocity. 

Using these expressions together with expression (5) tor the diffusion coef- 
ficient file transport equation (3) can be written 

2 • rp 


8n 

Tt 


„ f 6 2 n 2 +p 8n 2 * fp - 1) \ 

= kM " r + st + — "j 


+ 3.55 • 10' 6 (8000 * r -4 + 50) • 


((€’ •) - 
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(6) 
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The time evolution of the distribution function n can now be calculated by solv- 
ing Equation (6) with the appropriate initial - and boundary conditions. 

In order to obtain the initial distribution, the observed integral proton 
intensities have to be related to n. Nakada and Mead (1965) have shown that 

n (r, /z, J = 0, t) « r • j (r. E* *) < 7 > 

where j(r, E, t) is the differential flux (protons/cm 2 -s-sr-MeV). The inte- 
gral invariant J is equal Lo zero for particles mirroring at the equator. 

As shown previously by Davis and Williamson (1966) the integr al energy 
spectra measured by the scintillator experiment on board Explorer 26 can quite 
nicely be fitted with an exponential expression. Thus the differential flux is 
approximated by 

i( r o* E ’t 0 ) = " V E o * exp(~E/E 0 ) . (8) 

Using this equation and the definition of the magnetic moment of particles 
mirroring at the equator, 

fi = E/B (9) 

we obtain the following expression for n , 

n (V ^ J = °- t o) = r/E 0 exp(- M B/E 0 ) . (10) 

C is a constant of proportionality and the subscript on r and t indicate that n 
is evaluation at a radial distance r Q at a particular time t 0 . 
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The left-hand side of Figure 1 shows the measured radial intensity distribu- 
tion j(>E. r, t 0 ) of protons mirroring at the equator for the eight energies listed 
in the figure on day 111 (April 21) of 1965. From this mapping the distribution 
n (r. J - o, t 0 ) is obtained and shown on the right-hand side of the figure. 
Quadiaiic interpolation is used in the r -direction to obtain n in the region of 
fi-T space needed to solve Equation (6) numerically. Having obtained the initial 
distribution one has to choose the boundary conditions in order to solve Equation 
(6). Due to the increased atmospheric density and the corresponding heavy loss 
of particles from charge exchange and coulomb energy degradation, one knows 
and observes that the distribution function n goes to zero close to the earth. 
Somewhat arbitrarily n is set equal to zero at r - 1.5 earth radii. At r - 5.5 
earth radii the protea intensity is assumed to be independent of time at all 
energies, as shown by Soraas and Davis (1968) to be approximately true. They 
show that the intensities at L = 5.0 both before and after the April IS storm on 
the average run fairly constant, though significant variations occurred on a short 
time scale. For the time period and space region considered here, the results 
of the calculations do not depend critically on the boundary conditions as will be 
shown later. 

RESULTS AND DISCUSSIONS 

The time dependence of the protons have been calculated by solving Equation 
(6) numerically on a computer with the initial ai l boundary values as outlined 
in the previous section. 
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Three sets of values for the parameters in the diffusion coefficient expression 
(5) have been considered; 

Set A. k = 0.31 x 10” 9 , m = 0.0, and p = 10 as evaluated by Kakada and 
Mead (1965) using the observed frequency and size of sudden com- 
mencements and sudden impulses. 

Set B. m = 0.0 and p = 10 with k varied to minimize the RMS between cal- 
culated and measured values. 

Set C. All three parameters k, m, and p varied to minimize the RMS. 

With Set B and C where the least squares method was used to determine 
parameter values only data at L - 3.0, 3.5, and 4.0 earth radii were used. Data 
at L = 4.5 will be compared with calculated values but were excluded in the 
fitting process because the calculated values at this L-value are sensitive to the 
assumed L - 5.5 boundary values as is shown below. 

Figure 2 shows the calculated and observed integral intensities as a function 
of time for the four L-values 3.0, 3.5, 4.0, and 4.5 earth radii and for each of 
eight proton energies 134, 180, 220, 345, 513, 775, 1140, and 1700 keV. The 
calculated tim e variation for parameter set A are shown as dashed curves and 
for set B as solid curves. The measured points are the x's. The prestorm 
intensifies are indicated to show the rapid changes which took place during the 
storm and that the post-storm intensities slowly return to these values. It should 
be noted that the measured values and the values calculated from the initial dis- 
tribution on day 111 are not in complete agreement. This is because the meas- 
ured spectra are not exactly exponential as assumed. 


8 



Comparison between the calculated curves and measured points in Figure 2 
shows that the time variation computed with parameter set A does not follow the 
observed variation whereas those computed with parameter set B are reasonably 
good fits to the data. Thus the diffusion coefficient evaluated by Nakada and 
Mead (1965) is too small by a factor of 8 to account for the observed time varia- 
tion. This agrees with their further results that to get good agreement between 
the steady- state proton flux measured in 1962 and their calculated distribution, 
their diffusion coefficient needed to be increased by a factor of 8. Thus approxi- 
mately the same value of the diffusion coefficient is needed in order to bring the 
diffusion model in agreement with the steady state proton distribution measured 
in 1962 and with the time behavior of the protons after the April 18, 1965 geo- 
magnetic storm. This leads to the conclusion that the field fluctuations driving 
the diffusion had approximately the same intensity level during these two time 
periods. It also rather conclusively confirms that L-diffusion plays a major 
role in populating the outer zone protons. Of course other non- adiabatic proc- 
esses, such as that producing the rapid storm time changes observed by Soraas 
and Davis (1968,, are also important. The time variations calculated for set C 
are almost identical to the solid curves for set B in Figure 2. The diffusion 
coefficient parameters obtained for set B and C are listed in Table 1 along with 
the RMS of the fit measured in decibels. The units are such that D in Equation 
(5) is given in (earth radii) 2 per day when p. is in MeV/Gauss and r in earth 
radii. The standard deviations of the parameters listed are based on the 
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statistics of the tit. The value of m in set C came out slightly negative but not 
significantly different from zero. To determine how sensitive the RMS value is 
on the value of p the best fit to the experimental results was calculated, for a 
range of p values assuming n = 0. The RMS values exhibit a broad minimum 
for p values in the range 10 to 13. Due to the high correlation between the 
parameters k and p, possibly due to the limited and coarse coverage of data 
in the L -direction, it is not possible to determine p more accurately. The 
value of the diffusion coefficient D 2 in the L-range 3.0 to 4.0 is, however, 
fairly accurately determined. 

We notice that in Figure 2 the solid lines overlay the dashed lines for the 
three lowest energies at all L-values. This means that in this model the time 
behavior of these low energy protons are mainly governed by the loss mech- 
anisms and not by radial diffusion. 

The low energy protons at L-values above L = 3.5 are observed to decay 
faster during the time immediately following the storm time enhancement, than 
the model used here can account for. Losses due to charge exchange, which 
affect the low energy protons most, are expected to be of less importance at 
high L-values due to the reduction of the neutral hydrogen density with increas- 
ing radial distance from the earth. From this one should expect the decay rate 
of the low energy protons to go down with increasing L- value. It is thus believed 
that some other mechanism together with the ones considered here must control 
the decay rate of the low energy protons. Kennel and Petschek (1966) have 
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shown that pitch angle scattering by ion cyclotron noise sets an upper limit on 
protons which may be stably trapped. They have further shown that in the region 
about L = 4 the protons are near this limit- It is thus possible that the low’ 
energy protons exceeded this limit during their rapid storm -time enhancement. 
According to Kennel and Petschek (1966) the proton intensity should then be 
rapidly forced back to their stable trapping limit through a non-linear wave- 
particle interaction process. Tnis may then explain the fast decay of the low 
energy protons at high L- values following the storm-time enhancement. 

Another way of presenting the data is shown in Figure 3. The distribution 
function n is plotted vs. radial distance for different values of p the magnetic 
moment of the protons and for three different times, 0, 20 and 36 days after the 
storm. In the left-hand side of the figure n is plotted as obtained from the 
measured fluxes on Explorer 26. We can see how the shape of these curves 
changes as time progresses. The changes in n are more rapid at high 
L-values than at low. In this representation of the observed data the dif- 
fusive character of the storm time recovery becomes more apparent than 
in the plots of integral intensity vs. time. The right-hand side of the figure 
shows the calculated changes in n as computed from the diffusion equation 
using D 2 = 2.40 • 10" 9 r 10 (earth radii) 2 /day. ft is seen that the computed 
time evolution of the distribution function n in general follows the observed 
one quite closely. 
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To test how sensitive the diffusion calculations are on the boundary condi- 
tions at L = 5.5, Equation (6) was solved for different boundary conditions at 
this L-value using the diffusion coefficient D 2 = 2.4 * 10~ 9 r 10 (earth radii) 2 / 
day. The value of the distribution function n at L = 5.5 was reduced to half of 
the value used previously. In a time period of 45 days this change in the boundary 
conditions had no effects at L = 3.0 and 3.5. At L = 4.0 the calculation started 
to deviate after 20 days from the one presented in Figure 2. But still after a 
time of 45 days there was less than one db difference between the two calcula- 
tions. At L = 4.5 the two solutions started to deviate after 4 days and after 45 
days the difference was about 2 db. As an extreme the distribution function n 
was set equally to zero at L = 5.5. At L = 3.0, 3.5 and 4.0 the results of this 
calculation was essentially the same as in the case where n was reduced to half 
its value. At L = 4.5 though this last calculation showed a steady drop in the 
intensities after a 4 day period due to the presence of a sink one earth radii 
away. As the actual variations in the fluxes at 5.5 earth radii for most of the 
time are less than 50% and the variations take place on a time scale of a few 
days, the results of the diffusion calculations do not for the time period and 
space-region considered in this study, depend critically on that boundary 
condition. 

It is interesting to use the diffusion coefficient D 2 = 2.40 • 10" 9 r 10 (earth 
radii) 2 /day to calculate how the distribution function for the protons changes 
with time during a longer period using the same initial and boundary conditions 
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as before. In Figure 4 the result of this calculation is shown. The distribution 
function n is plotted vs. radial distance for two values of y, the magnetic 
moment of the protons. We can see how the shape of these curves changes as 
time progresses. The protons diffuse inwards approaching the steady state first 
at large radial distances and later closer to the earth. At first the changes are 
fairly rapid, but then they slow down as the spatial gradient in the distribution 

t 

function decreases. The time for the high energy protons which were depleted 
during the storm to complete 62% of the recovery to prestorm value are approx- 
imately 380 days at L = 3.5, 240 days at L = 4.0, 100 days at L = 4.5 and 14 
days at L = 5.0. The lower energy protons w'hich were enhanced reach steady 
state conditions faster due to their greater losses. In Figure 5 the calculated 
change in the spectral parameter E 0 for different L- values are shown vs. time. 
One can see how rapid E 0 changes and stabilizes at high L-values. The E 0 
values obtained after diffusion and loss processes have been working for 720 
days have the 1/L 3 dependence expected from a process conserving the two first 
adiabatic invariants and violating the third. 
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Table 1 


Set 

k • 10" 9 

m 

P 

RMS 

B 

(2.4 ± 0.1) 

0.0 

10 

0.86 

1 C 

0.73 ± 0.5 

-0.035 ± 0.7 

10.94 ± 0.2 

0.82 


Set B gives the value of the parameter k in Equation (5) giving the best fit to the 
data with m and p fixed. In set C the values of the parameters k, m, and p which 
give the best fit to the data allowing ail three to be varied are shown. The RMS 
of the fit is given in decibels. 
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FIGURE CAPTIONS 


Figure 1 - The initial distribution. The left-hand side of the figure shows the 
integral proton fluxes above various energies versus radial distance 
as measured after the April 18 storm on day 111 of 1965. The right- 
hand side of the figure shows the distribution function n for different 
values of the magnetic moment plotted vs. radial distance. 

Figure 2 - The time-behavior of the integral proton intensities at different L- 
values computed from the transport equation with two values of the 
diffusion coefficient, are compared with the experimentally measured 
values at L = 3.0, 3.5, 4.0 and 4.5. 

Figure 3 - The left-hand side of the figure shows the measured distribution 
function n, for constant magnetic moments plotted vs. radial 
distance for 0, 20 and 36 days after the storm. The right-hand 
side of the figure shows for the same days after the storm the 
distribution function n as calculated from the diffusion equation (6) 
sing D 2 = 2.40 • 10" 9 r 10 (earth radii) 2 /day plotted for constant 
magnetic moments vs. radial distance. 

Figure 4 - The time -evolution of the distribution function n calculated from the 
diffusion equation (6) using D 2 = 2.40 • 10“ 9 r 10 (earth radii) 2 /day 
plotted for constant magnetic moment vs. radial distance. The curves 
shown are for diffusion times of 0, 50, 150, 300 and 500 days. 



Figure 5 - The spectral e-folding energy E 0 at L-values 2.5, 3.0, 3.5, 4.0, 4.5 
and 1.0 computed from the solution of the transport equation and 
plotted vs. time. 
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Figure 1. The initial distribution. The left-hand side of the figure shows the integral proton fluxes above various energies versus 
radial distance as measured after the April 18 storm on day 111 of 1965. The right-hand side of the figure shows the distribution 
function n for different values of the magnetic moment plotted vs. radial distance. 
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Figure 2. The time-Lrohovior of the integral proton intensities at different L-values computed 
from the transport equation with two values of the diffusion coefficient, are compared with the 
experimentally measured values at L = 3.0, 3.5, 4.0 and 4.5. 
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Figure 3. The left-hand side of the figure shows the measured distribution function n, for constant magnetic moments plotted vs. 
radial distance for 0, 20 and 36 days after the storm. The right-hand side of the figure shows for the same days after the storm the dis- 
tribution function n as calculated from the diffusion equation (6) using Dj ~ 2.40 • 10” 9 r^® (earth radii)^/day plotted for constant mag- 
netic moments vs. radial distance. 
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Figure 4. The time-evolution of the distribution function n calculated from the dif- 
fusion equation (6) using Dj = 2.40 • 10“^ r*® (earth radii)^/day plotted for constant 
magnetic moment vs. radial distance. The curves shown are for diffusion times of 0, 
50, 150, 300 and 500 days. 
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